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Abstract. Extreme-mass-ratio inspirals (EMRIs), stellar-mass compact objects (SCOs) 
inspiralling into a massive black hole, are one of the main sources of gravitational waves 
expected for the Laser Interferometer Space Antenna (LISA). To extract the EMRI signals 
from the expected LISA data stream, which will also contain the instrumental noise as well 
as other signals, we need very accurate theoretical templates of the gravitational waves that 
they produce. In order to construct those templates we need to account for the gravitational 
backreaction, that is, how the gravitational field of the SCO affects its own trajectory. In general 
relativity, the backreaction can be described in terms of a local self-force, and the foundations 
to compute it have been laid recently. Due to its complexity, some parts of the calculation of 
the self-force have to be performed numerically. Here, we report on an ongoing effort towards 
the computation of the self-force based on time-domain multi-grid pseudospectral methods. 



1. Introduction 

Extreme-Mass-Ratio Inspirals (EMRIs) are one of the main source of gravitational waves for the 
planned Laser Interferometer Space Antenna (LISA) [1]. EMRIs are binary systems composed 
by a stellar compact object (SCO), like a neutron star or a stellar black hole, with masses 
in the range m = 1 — IO^Mq, inspiralling into a massive black hole (MBH), with masses in 
the range M = 10^ — 10 ''Mq, located at galactic centers. Then, the mass ratios involved are 
/u = m/M ~ 10"'^ — 10"^. During the inspiral phase the system is driven by the emission of 
gravitational radiation, and hence there is loss of energy and angular momentum that makes 
the orbit shrink until the SCO plunges into the MBH. It is expected that LISA will be able to 
detect 10 — 10'^ EMRI/yr up to distances with z < 1 [21 13] (see also |1]). In order to extract 
the signals produced by these systems, which will be buried in instrumental noise and the 
gravitational wave foreground (produced by compact binaries in the LISA band), and extract 
relevant physical information from them, we need to have an a priori theoretical knowledge of 
the gravitational waveforms with a high degree of precision. While techniques for constructing 
templates good enough for detection are getting ready, methods to build templates good enough 
for extraction of physical information are not yet fully developed. The main difficulty being 
that one requires a more precise treatment of the self-gravity of the SCO and its impact on the 
gravitational waveform. 

Frequency-domain approaches provide very accurate results in the computation of 
quasinormal modes and frequencies and also in the computation of the self-force for EMRIs 
with moderate eccentricities. However, the frequency domain approach has more difficulties 



when dealing with EMRIs with highly eccentric orbits, which are of interest for LISA, since 
one has to sum over a large number of modes to obtain a good accuracy, and convergence 
may be an issue. On the other hand, time-domain approaches are not affected much by the 
eccentricity of the orbit and may be more efficient for the case of high-eccentricity EMRIs. 
However, in the time domain numerical approach one has to deal with different physical scales 
(both spatial and temporal) present in the problem due to extreme mass ratios. Specifically, 
we have to computationally resolve large wavelength scales (comparable to the massive black 
hole) and also the scales in the vicinity of the SCO, which are crucial for evaluating the self- 
force. This means, in practice, that we have to resolve numerically large scales and small scales 
inside the same computational domain. Here we introduce a new time-domain scheme towards 
the computation of the self-force based on the PseudoSpectral Collocation (PSC) method (see, 
e.g. [5]). A remarkable property of this new scheme is that it avoids having to introduce in 
the problem a small spatial scale associated with the presence of the SCO (see, e.g. [6]). This 
is done by using multiple domains and locating the particle in the interface bewteen two of 
them. Therefore, we just need the resolution to describe the field near the particle, but not 
the particle itself, which makes the computation more efficient. Here, we describe how to apply 
this technique to the computation of the self-force on a charged scalar particle in circular orbits 
around a non-rotating black hole and show some results of the numerical implementation. 



2. Study of a Simplified Model 

We study a simplified model of an EMRI consisting of a particle with charge g, associated with a 
scalar field $(x^), orbiting a non-rotating MBH. In this model, the particle generates the scalar 
field, which in turn influences the particle trajectory (this effect can be described by a local force, 
called the self-force). Therefore, this model contains all the ingredients of the gravitational case 
(see details in [3 18] ) , where the particle motion is influenced by its own gravitational field instead 
of a scalar field. 

In the scalar field model the spacetime geometry is not dynamical. In our problem it has to 
describe the geometry of the MBH and therefore will be given by the Schwarzschild black hole 
metric: 

ds"^ = f{-dt^ + dr*'^) + r'^dn^ , dn^ = de"^ + sin^ edip'^ , (1) 

where (x^) = {t,r,6,(p) are the Schwarzschild coordinates, / = 1 — 2M/r, M is the MBH 
mass, and r* = r + 2Mln(r/2M — 1) is the tortoise coordinate. The scalar field equation is a 
wave-like equation with a singular source term due to the energy density of the particle's scalar 
charge (see, e.g. [9j) 

g'^'^VaVp^x) = -Attp = -Anq j dr 6^{x, z{t)) , (2) 

where x^ = z^{t) is the trajectory of the particle (without loss of generality, the orbital 
plane is taken to be the equatorial plane); denotes the associated canonical connection; 
r denotes proper time along the particle timelike worldline 7, and 5^{x, z{t)) is the invariant 
Dirac functional in Schwarzschild spacetime (see, e.g. P]) and describes the singular structure 
of the source term. The equations of motion for the particle, obtained from energy-momentum 
conservation, are 

m^=F'^ = g(g'^'^ + n'^n'^)V,$, = ^ , (3) 
or dr 

where m is the particle mass and its 4-velocity. The force diverges on the particle 
worldline, 7, due to the singularity that carries the particle field, and therefore it must be 
regularized |10] . An analysis of the solutions ofj2]) and ^ reveals (see [9] and references therein) 
that the field, <I>, can be split into two parts <I> = <I>^ -|- Where is the singular part. 



which contains the singular structure of the field and satisfies the same field equation, and ^> 
is the regular part, which satisfies an homogeneous wave equation: 
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The regular field is solely responsible of the deviation of the particle from geodesic motion around 
the MBH. Hence, the self-force which acts on the particle is = q{g^^ + u^^u'^)V ^^^^ . 

In order to solve the field equation Q we can make use of the spherical symmetry of the 
Schwartzchild metric and decompose ^ into scalar spherical harmonics. The equations for the 
different harmonic coefficients, r), decouple and satisfy 1+1 wave-like equations of the 

form 
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where V^(r') is the Regge- Wheeler potential for scalar fields on the Schwarzschild geometry and 
5™ is the harmonic coefficient of the singular source term generated by the particle, which is 
given by: 
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where an overbar denotes complex conjugation. The solution of (|5| is finite and continuous at 
the particle location, but it is not be differentiable in the sense that the evaluation of the radial 
derivative of the scalar field from the left and from the right of the particle yields different values, 
that is, there is a jump [see equation ( |12[ ) bellow]. Nevertheless, one can see that the contribution 
to the self-force of each harmonic is finite, it is the sum over i that becomes infinite. Taking 
advantage of this fact, a method to extract the singular part of the self- force mode by mode 
was developed, the so-called mode-sum regularization scheme |12 | I13 | [T^. Given the multipolar 
decomposition of the gradient of the scalar field 



'^^{x'^) = Y,'^iix'^), where ^iCx'^) = V„ ^ <^Tit,r)Y[ 
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each <I>^ is finite at the particle location. Using the mode-sum scheme, the regularized part of 
the gradient of the scalar field can be written as follows 
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and the multipoles of the singular part, , are known in a neighborhood of the particle 
worldline and are given by 
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where Aa, Ba, Ca, Da, ■■■ are called the regularization parameters [l3l HH |T5]. They are 
independent of i, but depend on the particle trajectory. The singular part corresponds to the 
first three terms which lead to quadratic, linear, and logarithmical divergences. The remaining 
terms form a convergent series that does not contribute to the self- force (each of them) . Since the 
only non-vanishing component of the regularization coefficients is the radial one, this is the only 
component of the self-force that is actually singular, and hence the only one to be regularized. 



3. Description of the Simulations 

We begin by defining the spatial domain where we want to obtain our solution: r* £ (—00, +cxd), 
where r* — > —00 corresponds to the horizon location (r = 2M), while r* — > +00 corresponds 
to spatial infinity. For computational reasons, we have to take a truncated radial domain (we 
could avoid this by compactifying the radial domain): r* G Q = [r* , r^], and hence we need to 
prescribe outgoing boundary conditions at the end points of this domain: 

d,{r ^Tit, rl)) - d,.. (r cl>™(t, rl)) = , a,(r r^)) + d,, (r cl>™(t, r^)) = . (10) 

On the other hand, in order to numerically compute the solutions, it is necessary to discretize 
our equations in time and space. We use the PSC for the spatial discretization, and the time 
discretization is performed by using a Runge-Kutta algorithm. Then, the harmonic components 
of the scalar field are approximated by an expansion in Chebyshev polynomials, 

N 

^?iv(i,r-*) = 5^a„(t)r„(X(r*)), Xe[-l,l], (11) 

n=0 

and the coefficients a„(t) are determined, in terms of a system of ordinary differential equations 
(ODEs), by forcing the equation ^ to be satisfied at a given set of collocation points 
{Xk}i^^Qi ^. For our calculations we took a Lobatto- Chebyshev grid, which includes the 
boundary points as nodes of the grid. The system of ODEs is solved by using a Runge-Kutta 
algorithm, typically a RK4 scheme. 

The singular source term ^ in our equation (jsj) would in principle spoil the exponential 
convergence of the PSC method. In order to avoid this, and also in order to avoid a small 
scale associated with the presence of the particle, we use a multi-domain scheme and locate 
the particle in the interface between two such domains. More specifically, we divide 0, into D 
sub-domains 17 = Lli[r*/",r*^] {i = 1, ...,D), with r*J\ = r*^, and evolve the equation (without 
the source term) in each sub-domain independently. In doing so we need to communicate the 
different subdomains, in particular the subdomains to the left and to the right of the particle, by 
using appropriate boundary conditions. These conditions can be derived from the equation ([5]), 
and can be written as junction conditions that describe the possible jumps in the solution 




In the case of circular orbits they are 



[^Tit,r*)], = 0, mTit,r*)]^ = 0, [d,.^Tit,r*)]^ = fSr. (12) 

Apart from these techniques, we also use the Fast Fourier Transform in order to evaluate 
derivatives of the scalar field in a quick way (see, e.g. [5]). Moreover, in order to reduce possible 
spurious high-frequency components of the numerical solution, we apply an exponential spectral 
filter to the solutions after every time step of the RK algorithm. 

4. Summary of Results 

We have developed a numerical code that implements the previous ideas in order to compute 
the self-force on a charged scalar particle in circular geodesies around a non-rotating MBH. 
In order to validite the code we have made a series of tests. In particular, we have evolved a 
Gaussian pulse in flat space (V; = 0) and also in the Schwarzschild geometry {Vi ^ 0), that is, 
we have evolved the equation ^ without a particle. We have done this using different numbers 
of collocation points, N, in order to check that we obtain the expected convergence of the PSC 
method for smooth solutions. In Figure [1] (top left) we show that indeed an estimation of the 
truncation error. In |a7v|, falls exponentially as the number of collocation points increases. 




Figure 1. Top left: Estimation of the truncation error, In |aAr| , versus the number of cohocation 
points, A^. The points he aroud a straight hne with negative slope, showing that the expected 
exponential convergence is obtained. The other figures show snapshots of the evolution of 
the (2, 2) mode (with rp = 6M) and the detailed structure near the particle. They show the 
continuity of {top right) and 9^$™ {bottom left), and the jump in dj.*^J^ {botton right). 



On the other hand, we have performed simulations including the charged particle in order 
to develop and test the communication between the subdomains, in particular between those 
around the particle where jumps in the solution can occur. In practice, this communication is 
implemented by using a penalty method, in such a way that the jump conditions in equation ( 12 ) 
are enforced in a dynamical way, by introducing new terms into the resulting ODEs that drive the 
solution towards the satisfaction of these conditions. These new terms have to be callibrated 
against the number of collocation points and the size of the time step. We have found that 
the penalty method, in combination with the use of a spectral filter, works very well and can 
produce accurate solutions to our equations. In order to illustrate this fact we have included 
three plots in Figure [T] corresponding to snapshots of the evolution of the (^, m) = (2,2) mode 
at the last stable circular orbit (r = 6M), which is probably the most demanding situation. At 
the top right, we show a snapshot of the scalar field, $2^) where it can be seen that the field 
it is continuous at the particle although it is not differentiable. At the bottom left we have a 
snapshot for the time derivative of the scalar field, d^^2^ which is also continuous (this only 
happens for the case of circular orbits [6]) but non-differentiable. Finally, at the bottom right, 
we show a snapshot of the radial derivative of the field, d^*^2i where we can see the ability 
our numerical techniques in resolving the jump accross the particle. In the small box inside this 
figure we zoom in near the particle location to show how well the method can deal with the 



sharp jump and at the same time produce smooth waves. 



5. Conclusions and Future Work 

We have presented a new time-domain computational framework to compute the self-force for 
EMRIs. We have implemented this new method numerically for the case of a charged scalar 
particular moving in circular geodesies around a non-rotating MBH. This method is based 
on the PSC method. One of the main ingredients of this technique is to consider a multi- 
domain grid and to locate the particle in the interface between two subdomains. In this way, 
we solve equations with smooth solutions inside each subdomain, which preserves the powerful 
convergence properties of the PSC method and, at the same time, we avoid introducing a small 
scale related to the size of the SCO. This in turn makes the algorithm quite efficient as compared 
to other treatments of the particle. The effects in the solution due to the presence of the particle 
source term appear in our method through the imposition of a set of jump conditions. Here, we 
have presented tests of the convergence of the numerical code we have built and we have also 
shown the ability of the code to resolve the field near the particle, including the jump in the 
radial derivative, and at the same time produce smooth waves. Calculations of the self-force on 
the particle are underway and will be presented elsewhere [IT] . 

In this work, we have restricted ourselves to the case of circular orbits, where the particle is 
always at the same radial position. In order to extend our calculations to generic eccentric orbits 
we need to introduce a method to deal with a moving particle, either by moving our domains or 
by changing the coordinate system so that the particle is always located at the same node. This 
is a project presently under development. Given that the astrophysically interesting EMRIs 
involve spinning MBHs, another important task for the future is to study how to implement all 
these techniques for the case in which the MBH is represented by the Kerr solution. 
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